classdef ArgoTempSaltProfileClass < handle
    % 用于储存Argo,XBT,GLD等温盐剖面数据的类
    % 一、构造函数使用方法
    % 1.通过构造函数构造一个空TempSaltProfileClass，然后通过setExtent设置范围
    % TSprofile = TempSaltProfileClass();
    % TSprofile.setExtent(extent);
    
    % 2.直接通过构造函数构建
    % TSprofile = TempSaltProfileClass('extent',extent);

    % 二、method介绍
    % 1.setExtent 同上 用于设置范围
    % 2.append 导入温盐剖面
    %%%% 实例
        % TSprofile.append(1, 1, 'sst', 1:10,'sss', 1:10, 'depth', 1:10, 'profileType','Argo','date',datetime(2011,1,1));
        % TSprofile.append(115, 15, 'sst', 1:10,'sss', 1:10, 'depth', 1:10, 'profileType','Argo','date',datetime(2011,1,1));


    properties
        % lon lat 
        % profile 保存sst sss depth
        % extent 设定范围
        % count 储存计数
        % ineddyArray 与eddy的位置关系
        lon
        lat
        profile
        extent
        ineddyArray
        count = 0;
    end
    
    methods

        function obj = ArgoTempSaltProfileClass(varargin)
            % 构造函数 初始化类 目前支持 extent初始化
            idx = find(strcmp(varargin,'extent'));
            if ~isempty(idx)
                obj.extent = varargin{idx + 1};
            end
        end

        function setExtent(obj, extent)
            obj.extent = extent;
        end

        function append(obj, lon, lat, varargin)
            % 如果存在位置限制
            if ~isempty(obj.extent)
                % [xbox, ybox] = area_box_set(obj.extent);
                % in = inpolygon(lon,lat,xbox,ybox);
                in = inpolygon(lon,lat,obj.extent(2,:),obj.extent(1,:));
                if in ~= 1
                    return;
                end
            end
                        
            idx_sst = find(strcmp(varargin,'temp'));
            idx_psss = find(strcmp(varargin,'sail'));
            idx_pres = find(strcmp(varargin,'pres'));
            idx_date = find(strcmp(varargin,'date'));
            idx_profileType = find(strcmp(varargin,'profileType'));
            if ~isempty(idx_psss) && ~isempty(idx_sst)
                % sst sss只要都不为空,视为输入了剖面
                obj.lon = [obj.lon; lon];
                obj.lat = [obj.lat; lat];
                obj.count = obj.count + 1;

                obj.profile(obj.count).pres = varargin{idx_pres+ 1};
                % obj.profile(obj.count).pres_adj = varargin{idx_pres+ 2};
                % obj.profile(obj.count).pres_qc = varargin{idx_pres+ 3};

                obj.profile(obj.count).date  = varargin{idx_date + 1};
                obj.profile(obj.count).profileType  = varargin{idx_profileType + 1};
                % 判断输入的剖面类型 如果idx为空则判断未输入
                if ~isempty(idx_sst)
                    obj.profile(obj.count).temp = varargin{idx_sst + 1};
                    % obj.profile(obj.count).sst_adj = varargin{idx_sst + 2};
                    obj.profile(obj.count).temp_qc = varargin{idx_sst + 2};
                end
                if ~isempty(idx_psss)
                    obj.profile(obj.count).sail = varargin{idx_psss + 1};
                    % obj.profile(obj.count).psss_adj = varargin{idx_psss + 2};
                    obj.profile(obj.count).sail_qc = varargin{idx_psss + 2};
                end
            else
                return;
            end
            function [xbox, ybox] = area_box_set(inside_extent)
                xboundary = [inside_extent(2,:)];
                yboundary = [inside_extent(1,:)];
                xbox = xboundary([1 1 2 2 1]);
                ybox = yboundary([1 2 2 1 1]);
            end

        end

        % 提取满足筛选要求的浮标例如有效深度和数据层数
        function [new_TS] = profileLimit(obj, varargin)
            idx_len = find(strcmp(varargin,'len'));
            idx_effictive_depth = find(strcmp(varargin,'effictive_depth'));
            idx_key = find(strcmp(varargin,'key'));
            selected_idx =[];
            for i = 1:length(obj.profile)
                depth = obj.profile(i).pres;
                obj.profile(i).effictive_Depth = depth(end);

                boolSelected = 1;
                if ~isempty(idx_len)
                    len_limit = varargin{idx_len + 1};
                    if length(depth) < len_limit
                       boolSelected = 0;                     
                    end
                end
                if ~isempty(idx_effictive_depth)
                    effictive_Depth_limit = varargin{idx_effictive_depth + 1};
                    if depth(end) < effictive_Depth_limit
                        boolSelected = 0; 
                    end
                end
                
                if ~isempty(idx_key)
                    key_limit = varargin{idx_key + 1};
                    if isnan(obj.profile(i).(key_limit))
                        boolSelected = 0; 
                    end
                end


                if boolSelected == 1
                    selected_idx = [selected_idx i];
                end
            end
            new_TS = ArgoTempSaltProfileClass.copy(obj, 'index', selected_idx);    
        end

        % 提取质量控制为好的浮标
        function [OkARgoTS, ErrorArgoTS] = error_qc_select(obj)
            idx_ok = [];
            idx_error = [];
            for i = 1:length(obj.profile)
                % pres_qc =  obj.profile(i).pres_qc;
                sst_qc = obj.profile(i).temp_qc;
                psss_qc = obj.profile(i).sail_qc;
                if all(psss_qc(:)==1)&&all(sst_qc(:)==1)
                    idx_ok = [idx_ok i];
                else
                    idx_error = [idx_error i];
                end
            end
            OkARgoTS = ArgoTempSaltProfileClass.copy(obj, 'index',idx_ok);
            ErrorArgoTS = ArgoTempSaltProfileClass.copy(obj, 'index',idx_error);
        end

         % ineddyArray中找到距离涡旋最近的XY
        function ineddydeltaXYDeal(obj)
            ineddy = obj.ineddyArray;

            for i = 1:length(ineddy)
                D = 9999;
                idx = 0;

                XY = ineddy(i).deltaXY;
                if isempty(XY)
                    continue;
                end
                [m,~] = size(XY);
                for j = 1:m
                    distance = sqrt(XY(j,1)^2 + XY(j,2)^2);
                    if distance < D
                        D = distance;
                        idx = j;
                    end
                end
                ineddy(i).deltaXY = XY(idx,:);
                type = ineddy(i).type;
                ineddy(i).type = type(idx);
                obj.ineddyArray = ineddy;
            end
        end

        % 对剖面进行各种物理量计算
        function profile_measure(obj)
            % for TS = 1:length(obj.lon)
            for TS = 1:length(obj.lon)
                TS
                TS_lon = obj.lon(TS);
                TS_lat = obj.lat(TS);
                value = fieldnames(obj.profile);
                %% 深度和压力检查
                if ~isfield(obj.profile,'depth')&&~isfield(obj.profile,'pres')
                    error('剖面缺少深度参数')
                elseif isfield(obj.profile,'depth')&&~isfield(obj.profile,'pres')
                    depth = obj.profile(TS).depth;
                    z = gsw_z_from_depth(depth);
                    obj.profile(TS).pres = gsw_p_from_z(z,TS_lat);
                else
                    pres = obj.profile(TS).pres;
                    z = gsw_z_from_p(pres,TS_lat);
                    obj.profile(TS).depth = gsw_depth_from_z(z);
                end
                
                %% 温度和盐度变换
                SP = obj.profile(TS).sail;
                T = obj.profile(TS).temp;
                p = obj.profile(TS).pres;
                p(p<0)=0;
                [SA, ~] = gsw_SA_from_SP(SP,p,TS_lon,TS_lat);

                
                CT = gsw_CT_from_t(SA,T,p);
                
                obj.profile(TS).pt = gsw_pt_from_t(SA,T,p,400);%p_ref = 400
                %% 位势密度和位势涡度
                rho = gsw_rho(SA,CT,p);
                f = gsw_f(TS_lat);
                drho = gradient(rho)./gradient(p);
                obj.profile(TS).rho = rho;
                obj.profile(TS).pv = f./rho.*drho;
                %% 计算动态高度
                try
                    obj.profile(TS).dh = gsw_geo_strf_dyn_height(SA,CT,p,1000);%p_ref = 1000
                catch
                    obj.profile(TS).dh = nan;
                end
            end
        end

        %移除计算结果
        function remove_measure_result(obj)
            try
                obj.profile = rmfield(obj.profile,'interp_sail');
                obj.profile = rmfield(obj.profile,'interp_rho');
                obj.profile = rmfield(obj.profile,'interp_pv');
                obj.profile = rmfield(obj.profile,'interp_dh');
                obj.profile = rmfield(obj.profile,'interp_pt');
    
                obj.profile = rmfield(obj.profile,'interp_pt_anomaly');
                obj.profile = rmfield(obj.profile,'interp_sail_anomaly');
                obj.profile = rmfield(obj.profile,'interp_rho_anomaly');
                obj.profile = rmfield(obj.profile,'interp_dh_anomaly');
                obj.profile = rmfield(obj.profile,'interp_pv_anomaly');
            catch
            end
        end
        % 批量进行插值方法
        function batchinterp(obj)
            [obj.profile] = ArgoTempSaltProfileClass.akimaInterpBatch(obj.profile, 'depth', 10:10:1000, 'sail');
            [obj.profile] = ArgoTempSaltProfileClass.akimaInterpBatch(obj.profile, 'depth', 10:10:1000, 'rho');
            [obj.profile] = ArgoTempSaltProfileClass.akimaInterpBatch(obj.profile, 'depth', 10:10:1000, 'pv');
            [obj.profile] = ArgoTempSaltProfileClass.akimaInterpBatch(obj.profile, 'depth', 10:10:1000, 'dh');
            [obj.profile] = ArgoTempSaltProfileClass.akimaInterpBatch(obj.profile, 'depth', 10:10:1000, 'pt');
        end


        % 通过减去背景场剖面方法求温盐异常
        function anomaly_measure_by_background(obj, backgroundArgo)
            obj.ineddydeltaXYDeal();
            InTwoR = arrayfun(@(x) x.inTwoR, obj.ineddyArray);
            InOneR = arrayfun(@(x) x.inOneR, obj.ineddyArray);
            In_idx = find(InTwoR>0 | InOneR>0);
            
            background_InTwoR = arrayfun(@(x) x.inTwoR, backgroundArgo.ineddyArray);
            background_InOneR = arrayfun(@(x) x.inOneR, backgroundArgo.ineddyArray);
            out_idx = find(background_InTwoR==0&background_InOneR==0);
            InArgo = ArgoTempSaltProfileClass.copy(obj,'index',In_idx);
            outArgo = ArgoTempSaltProfileClass.copy(backgroundArgo,'index',out_idx);

            out_date = arrayfun(@(x) x.date, outArgo.profile);
            for s = 1:length(InArgo.lon)
                t = InArgo.profile(s).date;
                out_idx = find(abs(days(out_date - t))<15);
                outArraylon = outArgo.lon(out_idx);
                outArraylat = outArgo.lat(out_idx);
                out_idxs = [];
                % out_ps = [];%添加权重
                for j = 1:length(out_idx)
                    d = sqrt((outArraylon(j) - InArgo.lon(s))^2 + (outArraylat(j) - InArgo.lat(s))^2)*111195;
                    if d<200000
                        out_idxs =  [out_idxs out_idx(j)];
                        % out_ps = [out_ps 1/(d)];
                    end
                end
                In_idx = s;
                anomly_measure(obj, outArgo,In_idx,out_idxs);
            end

            
            
            function anomly_measure(obj, outArgo,In_idx,out_idxs)
                if isempty(out_idxs)
                    obj.profile(In_idx).interp_pt_anomaly = nan;
                    obj.profile(In_idx).interp_sail_anomaly = nan;
                    obj.profile(In_idx).interp_rho_anomaly = nan;
                    obj.profile(In_idx).interp_dh_anomaly = nan;
                    return;
                end 
                % background_pt = 0;
                % background_sail = 0;
                % background_rho = 0;
                % background_dh = 0;
                % background_count = 0;
                % background_sum_p = 0;
                for s = 1:length(out_idxs)
                    background_pt(s,:) = outArgo.profile(out_idxs(s)).interp_pt;
                    background_sail(s,:) = outArgo.profile(out_idxs(s)).interp_sail;
                    background_rho(s,:) = outArgo.profile(out_idxs(s)).interp_rho;
                    if ~isnan(outArgo.profile(out_idxs(s)).interp_dh)
                        background_dh(s,:) = outArgo.profile(out_idxs(s)).interp_dh;
                    else
                        background_dh(s,:) = nan(size(outArgo.profile(out_idxs(s)).interp_pt));
                    end

                    % background_pt = background_pt + outArgo.profile(out_idxs(s)).interp_pt;
                    % background_sail = background_sail + outArgo.profile(out_idxs(s)).interp_sail;
                    % background_rho= background_rho + outArgo.profile(out_idxs(s)).interp_rho;
                    % background_dh= background_dh + outArgo.profile(out_idxs(s)).interp_dh;
                    % background_sum_p = background_sum_p + out_ps(s);
                    % background_count = background_count + 1;
                end
                [m,~] =size(background_pt);
                if m~=1
                    background_pt = nanmean(background_pt);
                    background_sail = nanmean(background_sail);
                    background_rho = nanmean(background_rho);
                    background_dh = nanmean(background_dh);
                end
                % background_pt =background_pt/background_count;
                % background_sail =background_sail/background_count;
                % background_rho =background_rho/background_count;
                % background_dh =background_dh/background_count;
                % background_pt =background_pt/background_sum_p;
                % background_sail =background_sail/background_sum_p;
                % background_rho =background_rho/background_sum_p;
                % background_dh =background_dh/background_sum_p;

                obj.profile(In_idx).interp_pt_anomaly = obj.profile(In_idx).interp_pt - background_pt;
                obj.profile(In_idx).interp_sail_anomaly = obj.profile(In_idx).interp_sail - background_sail;
                obj.profile(In_idx).interp_rho_anomaly = obj.profile(In_idx).interp_rho - background_rho;
                obj.profile(In_idx).interp_dh_anomaly = obj.profile(In_idx).interp_dh - background_dh;
            end
        end
        
        % 通过气候态数据求温度异常
        function anomaly_measure_by_car09(obj)
            for TS = 1:length(obj.lon)
                TS
                TS_lon = obj.lon(TS);
                TS_lat = obj.lat(TS);
                date = datetime(obj.profile(TS).date.Year, obj.profile(TS).date.Month, obj.profile(TS).date.Day);

                t = days(date - datetime(date.Year,1,1)) + 1;
                cars09file = ['H:\ftp\CARS09\数据处理\SCS_Cars09_Potential_temp\' num2str(t) '.mat'];
                % load(cars09file,'car_lon','car_lat','car_temp','car_sail','car_depth');
                load(cars09file,'car_lon','car_lat','car_pt','car_sail','car_rho','car_dh');
                TS_lon = double(roundn(TS_lon,-1));
                TS_lat = double(roundn(TS_lat,-1));
                lonI = abs(car_lon - TS_lon) < 0.0001;
                latI = abs(car_lat - TS_lat) < 0.0001;
                pt = squeeze(car_pt(lonI,latI,:));
                sail = squeeze(car_sail(lonI,latI,:));
                rho = squeeze(car_rho(lonI,latI,:));
                dh = squeeze(car_dh(lonI,latI,:)); 
                % pv = squeeze(car_pv(lonI,latI,:)); 
                % car_depth = car_depth';
                % %% 温度和盐度变换
                % SP = sail;
                % T = temp;
                % p = gsw_p_from_z(-car_depth,TS_lat);
                % [SA, ~] = gsw_SA_from_SP(SP,p,TS_lon,TS_lat);
                % CT = gsw_CT_from_t(SA,T,p);               
                % pt = gsw_pt_from_t(SA,T,p,400);%p_ref = 400
                % %% 位势密度和位势涡度
                % rho = gsw_rho(SA,CT,p);
                % %% 计算动态高度
                % dh = gsw_geo_strf_dyn_height(SA,CT,p,1000);%p_ref = 1000
                obj.profile(TS).interp_pt_anomaly = obj.profile(TS).interp_pt - pt';
                obj.profile(TS).interp_sail_anomaly = obj.profile(TS).interp_sail - sail';
                obj.profile(TS).interp_rho_anomaly = obj.profile(TS).interp_rho - rho';
                obj.profile(TS).interp_dh_anomaly = obj.profile(TS).interp_dh - dh';

                % %%% 通过rho计算pv
                % obj.pv_anomaly_measure();
            end
        end

        % function pv_anomaly_measure(obj)
        %     for TS = 1:length(obj.lon)
        %         rho = obj.profile(TS).interp_rho_anomaly;
        %         %%% 通过rho计算pv
        %         drhoz = gradient(rho);
        %         dz = gradient(obj.profile(TS).interp_depth);
        %         drho = drhoz./dz;
        %         f = gsw_f(15);
        %         pv = f./rho.*drho;
        %         obj.profile(TS).interp_pv_anomaly = pv;
        %     end
        % end




        % 检查错误剖面 目前主要检查错误的位温异常剖面
        %一般没用特殊情况使用
        function [new_TS] = anomly_check(obj)
            selectedindex = [];
            for i = 1:length(obj.ineddyArray)
                boolselected = true;
                XY = obj.ineddyArray(i).deltaXY;
                type = obj.ineddyArray(i).type;
                X = XY(1);
                Y = XY(2);
                d = sqrt(X^2 + Y^2);
                %首先不能有nan值
                if isnan(obj.profile(i).interp_pt_anomaly)
                    warning(['已剔除序号为第' num2str(i) '个剖面'])
                    boolselected=false;
                end
                %其次不能出现1.1半径内出现大面积负异常情况
                if d <=1.1 & ~isnan(obj.profile(i).interp_pt_anomaly)
                    interp_pt_anomaly =  obj.profile(i).interp_pt_anomaly;
                    if interp_pt_anomaly(30) * type > 0 
                        warning(['已剔除序号为第' num2str(i) '个剖面'])
                        boolselected=false;
                    end
                end
                %最后对高深度位温异常不能超过1°
                if abs(obj.profile(i).interp_pt_anomaly(80))>1
                        warning(['已剔除序号为第' num2str(i) '个剖面'])
                        boolselected=false;
                end

                %低深度位温异常不能超过5°
                if abs(obj.profile(i).interp_pt_anomaly(10))>5
                        warning(['已剔除序号为第' num2str(i) '个剖面'])
                        boolselected=false;
                end

                if boolselected==true
                    selectedindex = [selectedindex i];
                end
            end
            new_TS = ArgoTempSaltProfileClass.copy(obj, 'index', selectedindex);    
        end
        % 消除剖面系统误差
        function [error] = remove_system_error_by_key(obj,key)
            TSMEC = TSprofileMatchEddyClass(obj,[]);
            [Untrapped_Argo] = TSMEC.profileDetectBymultiR("outTwoR");
            Untrapped_Argo.data_clean_by_key(key);
            [error] = system_error_measure(Untrapped_Argo, key);
            remove_error(obj,error,key);
            function [error] = system_error_measure(Untrapped_Argo, key)
                m = length(Untrapped_Argo.lon);
                n = length(Untrapped_Argo.profile(1).(key));
                for Untrappedi = 1:m
                    Matrix(:,Untrappedi) = Untrapped_Argo.profile(Untrappedi).(key);
                end
                for Untrappedj = 1:n
                    error(Untrappedj) = nanmean(Matrix(Untrappedj,:));
                end
            end
            function remove_error(obj,error,key)
                m = length(obj.lon);
                n = length(obj.profile(1).(key));
                for Ri = 1:m
                    Matrix(:,Ri) = obj.profile(Ri).(key);
                end
                for Rj = 1:n
                    Matrix(Rj,:) = Matrix(Rj,:) - error(Rj);
                end
                for Ri = 1:m
                    obj.profile(Ri).(key) = Matrix(:,Ri);
                end

            end

        end

        %四分位法数据清洗
        function data_clean_by_key(obj,key)
            m = length(obj.lon);
            n = length(obj.profile(1).(key));
            for i = 1:m
                Matrix(:,i) = obj.profile(i).(key);
            end
            for j = 1:n
                Matrix(j,:) = data_clean(Matrix(j,:),'std');
            end
            
            for i = 1:m
                obj.profile(i).(key) = Matrix(:,i);
            end
            
            function data = data_clean(data,method)
                if strcmp(method,'quantile')
                    p = [0.25 0.5 0.75];%分位数位置向量
                    q = quantile(data,p);
                    IQR = q(3) - q(1);%计算四分位距
                    lower_bound = q(1) - 1.5*IQR;%下界
                    upper_bound = q(3) + 1.5*IQR;%上界
                    data(data<lower_bound|data>upper_bound) = nan;%出界数据舍弃；
                else 
                    if strcmp(method,'std')
                        data = Three_std_data_clean(data);
                        data = Three_std_data_clean(data);
                        data = Three_std_data_clean(data);
                    end
                end
                
                function data = Three_std_data_clean(data)
                    mean = nanmean(data);
                    std = nanstd(data);
                    data(data<mean - 3*std|data>mean + 3*std) = nan;
                end
            end
        end

    end
    methods(Static)

        function [profile] = akimaInterpBatch(profile, interpZkey, interpZ, interpdatDatakey)
            for i = 1:length(profile)
                x = profile(i).(interpZkey);
                y = profile(i).(interpdatDatakey);
                try
                    profile(i).(['interp_' interpZkey]) = interpZ;
                    interpResult = makima(x,y,interpZ);
                    profile(i).(['interp_' interpdatDatakey]) = interpResult;
                    profile(i).effictive_Depth = x(end);
                catch
                    disp(num2str(i));
                    profile(i).(['interp_' interpdatDatakey]) = nan;
                    profile(i).interp_Depth = nan;
                    profile(i).effictive_Depth = nan;
                end
            end
        end
        
        % 需要插值才能使用
        function [xi,yi,zi,vq] = scatterInterp(profile, XY, interpXY, key,varargin)
%%%%%%%%%%%%%%%%%
            zi = 10:10:1000;
            xi = interpXY(1,:);
            yi = interpXY(2,:);
            [xi,yi,zi] = ndgrid(xi,yi,zi);
            x = [];
            y = [];
            z = [];
            f = [];
            for i = 1:length(profile(1).(key))
                i
                x_item = double(XY(1,:)');
                y_item = double(XY(2,:)');
                z_item = ones(length(x_item),1)*zi(1,1,i);
                f_item = double(arrayfun(@(x) x.(key)(i), profile)');

                TF_idx = find(~isnan(f_item));
                x_item = x_item(TF_idx);
                y_item = y_item(TF_idx);
                z_item = z_item(TF_idx);
                f_item = f_item(TF_idx);

                x = [x;x_item];
                y = [y;y_item];
                z = [z;z_item];
                f = [f;f_item];
            end
            mask = ones(size(xi));
            pm = ones(size(xi)) / (xi(2,1,1)-xi(1,1,1));
            pn = ones(size(xi)) / (yi(1,2,1)-yi(1,1,1));
            po = ones(size(xi)) / (zi(1,1,2)-zi(1,1,1));
            idx_len = find(strcmp(varargin,'len'));
            if ~isempty(idx_len)
                len = varargin{idx_len + 1};
            else
                len = 0.1;
            end
            idx_lambda = find(strcmp(varargin,'lambda'));
            if ~isempty(idx_lambda)
                lambda = varargin{idx_lambda + 1};
            else
                lambda = 50;
            end
            vq = divand(mask,{pm,pn,po},{xi,yi,zi},{x,y,z},f,len,lambda);
%%%%%%%%%%%%%%%%%%%
            % for i = 1:length(profile(1).(key))
            %     % i = 20
            %     i
            %     x = double(XY(1,:)');
            %     y = double(XY(2,:)');
            %     f = double(arrayfun(@(x) x.(key)(i), profile)');
            %     % f = arrayfun(@(x) x.(key)(i), profile, 'UniformOutput', false);
            %     % for i = 1:length(f)
            %     %     v(i) = f{i};
            %     % end
            %     % v = double(v');
            %     % f = v;
            %     TF_idx = find(isoutlier(f) == 0 | isnan(f));
            %     x = x(TF_idx);
            %     y = y(TF_idx);
            %     f = f(TF_idx);
            % 
            %     [xi,yi] = ndgrid(interpXY(1,:), interpXY(2,:));
            %     mask = ones(size(xi));
            %     pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
            %     pn = ones(size(xi)) / (yi(1,2)-yi(1,1));
            %     % fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,15,100,'factorize',0);
            %     % fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,5,20,'factorize',0);
            %     fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,2.2,20,'factorize',0);
            %     vq(:,:,i) = fi;
            % end    
        end

        function [xx,yy,zz,vq] = grid_scatterInterp(profile, XY, interpXY, key,varargin)
            zi = 0:10:1000;
            xi = interpXY(1,:);
            yi = interpXY(2,:);

            x = [];
            y = [];
            z = [];
            f = [];
            for i = 1:length(profile(1).(key))
                i
                x_item = double(XY(1,:)');
                y_item = double(XY(2,:)');
                z_item = ones(length(x_item),1)*zi(i);
                f_item = double(arrayfun(@(x) x.(key)(i), profile)');

                TF_idx = find(~isnan(f_item));
                x_item = x_item(TF_idx);
                y_item = y_item(TF_idx);
                z_item = z_item(TF_idx);
                f_item = f_item(TF_idx);

                x = [x;x_item];
                y = [y;y_item];
                z = [z;z_item];
                f = [f;f_item];
            end
            F = scatteredInterpolant(x,y,z,f);

            [xx,yy,zz] = meshgrid(xi,yi,zi);
            vq = F(xx,yy,zz);
        end

        function [xx,yy,zz,vq] = cressman_scatterInterp(profile, XY, interpXY, key,varargin)
            xi = interpXY(1,:);
            yi = interpXY(2,:);
            zi = 10:10:1000;
            [xx,yy,zz] = meshgrid(xi,yi,zi);
            for i = 1:length(profile(1).(key))
                % i = 20
                i
                x = double(XY(1,:)');
                y = double(XY(2,:)');
                f = double(arrayfun(@(x) x.(key)(i), profile)');
                [vq(:,:,i)] = cressman_interp(xi,yi,x,y,f,0.1,2);
            end   

            function [DATA] = cressman_interp(lon,lat,data_lon,data_lat,data,R,l_w)   % in degree of longitude
            % lon : new 1D longitude grid.
            % lat : new 1D latitude grid.
            % l_w=input('Which weighting function? (1=Square, 2=Exponential) ');
            % init_b : we need an initial backgroud data in the grid to be interpolated.
            
                init_b=zeros(length(lon),length(lat));
                Rsq=R*R;
                anal=zeros(length(lon),length(lat));
                if l_w == 1
            %% Square weighting function
                    for ii=1:length(lon)
                        for jj=1:length(lat)
                            total_weight=0;
                            for kk=1:length(data)
                                dx=abs(lon(ii)-data_lon(kk));
                                dy=abs(lat(jj)-data_lat(kk));
                                dis=sqrt(dx.^2+dy.^2);
                                if (dis<R)
                                    weight = (R.^2-dis.^2) ./ (R.^2+dis.^2);
                                    anal(ii,jj) = anal(ii,jj) + (data(kk)-init_b(ii,jj))*weight;
                                else
                                    weight=0;
                                end
                                total_weight=total_weight+weight;
                            end
                            if total_weight == 0
                                anal(ii,jj)=0;
                            else
                                anal(ii,jj)=anal(ii,jj)./total_weight+init_b(ii,jj);
                            end
                        end            
                    end
                    DATA=anal;
            %% Exponential weighting function
                else    
                   for ii=1:length(lon)
                        for jj=1:length(lat)
                            total_weight=0;
                            for kk=1:length(data)
                                dx=abs(lon(ii)-data_lon(kk));
                                dy=abs(lat(jj)-data_lat(kk));
                                dis=sqrt(dx.^2+dy.^2);
                                if (dis<R)
                                    weight = exp(-dis.^2./ (2*R.^2));
                                    anal(ii,jj) = anal(ii,jj) + (data(kk)-init_b(ii,jj))*weight;
                                else
                                    weight=0;
                                end
                                total_weight=total_weight+weight;
                            end
                            if total_weight == 0
                                anal(ii,jj)=0;
                            else
                                anal(ii,jj)=anal(ii,jj)./total_weight+init_b(ii,jj);
                            end
                        end            
                    end
                    DATA=anal;
                end
            end
        end

        function obj = copy(oldobj, varargin)
            obj = ArgoTempSaltProfileClass();
            TS_Index_idx = find(strcmp(varargin,'index')); 
            if isempty(TS_Index_idx)
                TS_Index = 1:length(oldobj.lon);
            else
                TS_Index = varargin{TS_Index_idx + 1};
            end
            obj.lon = oldobj.lon(TS_Index);
            obj.lat = oldobj.lat(TS_Index);
            obj.profile = oldobj.profile(TS_Index);
            if ~isempty(oldobj.ineddyArray)
                obj.ineddyArray = oldobj.ineddyArray(TS_Index);
            end
            obj.count = length(TS_Index);
            obj.extent = oldobj.extent;
        end
    end
end

